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Q ■ Abstract 



We introduce economic models based on Boolean Delay Equations: this 
formalism makes easier to take into account the complexity of the interactions 
between firms and is particularly appropriate for studying the propagation of 
an initial damage due to a catastrophe. Here we concentrate on simple cases, 
which allow to understand the effects of multiple concurrent production paths 
as well as the presence of stochasticity in the path time lengths or in the network 
structure. 

In absence of flexibility, the shortening of production of a single firm in an 
isolated network with multiple connections usually ends up by attaining a finite 
fraction of the firms or the whole economy, whereas the interactions with the 
outside allow a partial recovering of the activity, giving rise to periodic solu- 
tions with waves of damage which propagate across the structure. The damage 
propagation speed is strongly dependent upon the topology. The existence 
of multiple concurrent production paths does not necessarily imply a slowing 
down of the propagation, which can be as fast as the shortest path. 



1 



1 Introduction 



In most economic models, the production system is modeled as a unique representative 
producer (e.g., with the Cobb-Douglas function of the Solow model) or as a set of 
sectors, in which there is a unique representative producer by sector (e.g., in General 
Equilibrium Models). But the real production system can best be seen as a network 
composed of firms, producing different goods and services, and connected by suppliers 
and customer links. In such a network, the firm j supplies a fraction of its production 
to other firms, that use this production as an input for their own production function 
(see [Fig. 1], for instance). 

Introducing the role of networks in the economic system can lead to complex 
endogenous economic dynamics [Tj. But the network formalism is also well adapted 
to study the cascade effects generated by exogenous events, positive in the case of 
new orders from the market (2j [3j HJ El El U\ IE] and negative in the case of financial 
crisis [9], random local strikes affecting production [HUE], or natural disasters [121 
E1E]- I n classical economic models, where the production system is represented as 
a unique representative producer or as a small set of representative producers, indeed, 
exogenous events effects on the numerous firms have to be averaged over all firms, or 
over all firms of each sector. Because such effects are often highly heterogeneous, and 
because consequences and responses are highly nonlinear, this averaging process can 
bias the analysis, and makes it impossible to assess the consequences of an exogenous 
shock. 

The case of natural disasters is particularly interesting, because disaster impacts 
are very heterogeneous and affect especially strongly a small set of firms [T5| fT6]. In 
such cases, the total economic impact of the catastrophe can be much higher than 
the direct impacts of the event, because indirect effects through supply chains can be 
large. For instance, an earthquake that destroys a bridge can cause losses that are 
much larger than the value of the bridge, because impacts on transportation costs 
and duration can impair production in many firms. 

Such results have been reported using input-output models at the sector level 
[T71 [T8l [T9l [T3j . or specific network models [20]. To account for heterogeneity at the 
firm level, an appropriate input-output formalism has been developed in [121 El E] , 
and used to analyze disaster consequences. This approach has demonstrated that the 
shape and the structure of the network play an important role in disaster vulnerability, 
justifying the introduction of network effects in the economic assessment of natural 
disasters. The purpose of the present manuscript is to simplify as much as possible 
the formalism, within the framework of Boolean Delay Equations (BDEs), to be able 
to go beyond the simple observation that the network structure has an influence on 
disaster consequences, and to analyze this dependency. 

BDEs are semi-discrete dynamical systems, whose discrete variables evolve in 
continuous time; they have been introduced about 25 years ago [2U 1221 123], and they 
are related to the kinetic logic of Thomas [21] : nevertheless, in BDEs, the memory of 
the system can contain more and more information as the time goes on, allowing for 
solutions of increasing complexity, which hence display deterministic chaotic behavior, 
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as it has been recently observed experimentally [25]. Apart from their intriguing 
mathematical properties, BDEs represent an useful tool in the modeling of complex 
systems, characterized by threshold behavior, multiple feedbacks, and distinct time 
delays. They have been in particular successfully applied to the study of climate 
dynamics [26], of earthquake physics [271 EEL an d of real life problems [291 EH] (see 
[3Tj for a recent review). The present work is a first step towards the application of 
BDEs to phenomena in the economic realm and, at the same time, it deals with the 
clarification of the role played by stochasticity in these systems. 

By using BDEs, we plan to take into account the importance of the topology of the 
network, and of the different delays involved in the production paths [21 El SI El El Ej , 
on the total losses due to the initial catastrophe. The manuscript will develop the 
role of the connectivity of the structure: how the multiplicity of products needed by 
each industry affects the vulnerability of the whole industry, and what are the effects 
of multiple production paths from one firm to down-stream production firms. We are 
specially interested by long term production shortage spatio-temporal patterns gener- 
ated by the generic asynchrony due to variable time lengths of concurrent production 
paths. 

Even though we are using Boolean variables to describe the production of indi- 
vidual firm, a zero value should not be interpreted in terms of a destruction of the 
production unity: we rather mean that some shortage has been generated through 
production interactions; similarly, a level of one simply implies that the firm has 
recovered from previous impairing of production. When taking instead the interpre- 
tation of the zero value as a complete destruction of the firm, our models could be 
interesting from the point of view of the study of damage spreading in an Opera- 
tional Research framework. Finally, we will show that a macroeconomic quantity, 
measuring the density of fully active firms, allows to characterize the propagation of 
the consequences of the initial catastrophe in the various considered cases. 

2 The models 

A realistic representation of the economy at the firm level would make possible to 
assess both the direct and the indirect losses due to a catastrophic event, by taking 
into account also the backward and forward propagations of the event's consequences, 
hence the failure avalanches and the ripple effects which move across the chains of 
suppliers and producers in the network. For instance, the analysis of the impacts on 
the regional economy of the Northridge earthquake [T5l 132] makes it evident that a 
catastrophe can have very heterogeneous repercussions, and that the indirect losses 
due in particular to damages to the transport infrastructure system can be definitely 
higher than the direct losses themselves. Similar conclusions are obtained by the 
study of the consequences of the Loma Prieta and Northridge earthquakes [16], which 
underlines the importance of taking into account both the direct and the indirect 
losses, too; moreover this work shows that the repercussions can be less extended on 
firms belonging to a market larger than the strictly local one, hence suggesting to 
consider also economic models with adaptation (see, e.g., [20]). 
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Whereas the first steps towards a model suitable for realistic evaluations have 
been presented in [121 E2 E3], we introduce here some quite simplified toy models, 
which we will study within the framework of Boolean Delay Equations: the present 
approximation neglects some key ingredients, such as in particular the division in 
sectors of the economy, and the related observation that a given firm can usually 
choose between more providers of the same good, also in a strictly local economy 
Nevertheless, with these over-simplifications, the numerically observed behaviors can 
be understood in detail from the theoretical point of view, and we find that our results 
are a good starting point for future analysis. 

We take the Boolean variable Xi(t) = to mean that firm % at time t is impaired 
and can not fully produce (respectively, Xi(t) = 1 implies that it is not impaired): 
this can be due both to the firm being itself damaged or to the fact that it lacks 
of the necessary inputs, since some of its suppliers (or some of the suppliers of its 
suppliers and so on) were previously damaged. Therefore, this simple modeling takes 
into account the role of the chains of suppliers and producers in real economies. 

We study networks of N firms, assumed to be placed on the vertexes of a directed 
graph defined by the connectivity matrix A. In our notation, = 1 if part of 
the output of firm j (the supplier) is needed as input for firm i (the customer) and 
Aij = otherwise. As a first step for assessing the losses due to the propagation of the 
consequences of a catastrophic event, hence the vulnerability of an economic network 
to a natural disaster, we analyze the vulnerability of the present models to the initial 
damage of a single firm. Therefore we will assume, without loss of generality, that the 
firm in the origin is initially destroyed during an interval of time of length r c . We will 
look in particular at the dependence of the propagation of this event on the structure 
of the matrix A, hence on the resulting network topology, and at its dependence on 
the ensemble of the delays that exist in the economic system. 

We consider both isolated networks, the free models, and networks interacting 
with the outside, hence economies with adaptation, the forced models. In the first 
cases, the availability Sji of the good manufactured by firm j at firm i is simply 
delayed from production Xj with a constant time length according to: 

Sji(t) = Xj(t — Tij) free models. (1) 

In the presence of adaptation, which means that the economy is not locally isolated 
but there is instead the availability of some external rescue input, the stock is 
made again disposable after that the production of the firm i has been impaired for 
a duration of time that we take, in a first approximation, still equal to r^j: 

Sjiit) =Xi(t — Tij) V Xj(t — Tij) forced models, (2) 

where V is the OR operator and (■) means the Boolean negation. These equations 
give the truth table reported in [Tab. 1]. 

In our simple toy models, the production Xi(t), of firm i at time t, needs the 
presence of all the stocks of goods usually provided by the suppliers {j : A^ ^ 0} 
to which the firm i is connected in the network. Besides, we assume that Xi(t) also 
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j and i inactive, the stock is supplied from outside 
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j active and i inactive, the good is stocked 
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j inactive and i active, the stock is finished 


1 


1 


1 


j active and i active, the stock is updated 



Table f : The output table of the stock Sji of a given product as Boolean function of 
the activities Xj, Xj of the customer firm i and of the supplier firm j, in the free models 
described by Eq. ([I]) and in the forced models described by Eq. ([2]), respectively. 

depends upon the value of an external coefficient, which is an independent 

Boolean variable that allows in particular for the initial destruction of the firm in the 
origin. These constraints can be expressed, both for the free and the forced models, 
by the system of N BDEs: 

N 

Xi(t) = m(t) Y[A i3 V Sji(t) for i = 1, . . . , N. (3) 
j'=i 

Here the product fj, which runs over all the sites of the network, means the Boolean 
operator AND, A. The dependence on the delays {nj} is implicit, through the stocks, 
given by Eq. ([!]) or Eq. (j2J). 

In order to get well defined solutions for this system, one has to fix the initial 
values of the set of variables {xi(t)} in the interval [— r max , 0], where r max is the largest 
possible delay. Moreover, one has to choose the behaviors of the external functions 
{(ii(t)}. We only study the simple cases in which the economy starts undamaged 
and all the external functions take the value one during the evolution, apart from the 
destruction of the single firm in the origin, i.e. on the node i — 1. Notice that this 
firm is not definitely eliminated form the network: we assume instead that it is forced 
to stop the activity from the time t = up to the time r c ; this means H\{t) = for 
te [0,r c ]. 

We compare results for deterministic delays, chosen all equal to tq, with those 
for the more realistic situation of random delays {Tij}, independently and uniformly 
distributed in the interval [r min , T max \, with: 

V(Tij) = 5 (nj/Tmin - k) , (4) 

k=l 

i.e. the random variables are multiples of the unity of time fixed by r m j n . We usually 
take To = r min = 1 day and r max = 10 days. Though the possibility of defining BDE 
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Fr, V(n d ) 
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Fo, V{r id ) 


CM, n = 1 
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CM, n > 1 


3.4 


3.5 


3.7 


3.8 


DRG 


4.2 


4.3 


4.4 


4.4 



Table 2: We summarize here the main characteristics of the models studied in the 
following, and we indicate in each box the section where the model with all the 
corresponding details is discussed and the figures showing our results are presented. 
We label Fr the free models described by Eq. ([1]) and Fo the forced models obeying 
to Eq. (j2J); then, we look both at the case in which all the delays are deterministically 
taken equal to the same time unit To, and at the random delay one, with Vfcj) 
defined by Eq. (j4]); still, CM refers to the deterministic braid chain network obtained 
from a circulant matrix, and DRG to the directed random graph, with V(Aij) given 
by Eq. (jSJ; finally, we distinguish between single (n = 1), and multiple (n > 1), 
input-output connections in the braid chain, whereas the models defined on the DRG 
with average in/out-degree z both smaller and larger then one are treated in the same 
sections. 



systems characterized by random delays with a given probability distribution was 
already stated in [2T], the present study is the first example in this direction to our 
knowledge. 

Still, in the first part of this work we look at a deterministic braid chain network, 
corresponding to a connectivity matrix A with circulant structure (see [Fig. 1]), 
whereas in the second part we study a directed random graph (DRG). We will consider 
in particular the family of directed random graphs, D(N,p) [331 131], obtained by 
generalizing the rule of the well known Erdos-Renyi undirected random graph (RG) 
[33| I3SJ EZ], where each of the N(N — 1) directed links is present or absent with 
the same independent probability p £ [0, 1]. Hence, the elements of the connectivity 
matrix {A^} are random variables, assumed to be independently and identically 
distributed: 

V(A \ - S S ( Al ^ for i = J 

^ \ p6 (A i:j - 1) + (1 - p)S (A^) for i ^ j 

The probability p is straightaway related to the mean number of input-output con- 
nections (i.e. the in/out- degree), z = (k in ) = (k out ) = (N — l)p ~ Np, of the resulting 
directed network. 

In the considered context of BDEs with stochasticity, the delays {t^-} and the 
elements of the connectivity matrix {A^} are quenched random variables: when 
defining the BDE system one fixes their values, with probabilities given by Eq. (j3J) 
and Eq. (jSJ), to a given disorder configuration, and they are assumed to be constant 
on the time scale of the evolution we are interested in. 

With the aim of describing the solutions, and particularly their behaviors in the 
limit of a large number N of variables, we find interesting to introduce a macroeco- 
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nomic observable, that is to say the average density of firms which maintain their 
production, pit). In the deterministic case it is given by: 



1 N 



In the presence of stochasticity, in the delays and/or in the elements of the connectiv- 
ity matrix, the behaviors of the systems are better captured from the density averaged 
over disorder: 

(p(t)) = I \{dr^V{r id ) I Y[dA hk V(A hk )p(t), (6) 

i,j h,k 

that can be computed numerically by considering an enough large number, J\f s , of 
different disorder configurations. We also look at the density averaged over a time 
window of length corresponding to the smallest possible delay, i.e. r — r in the 
deterministic case and r = T rnin in the stochastic one: 

/t+T 
p(s)ds. (7) 

Finally, in the considered simple toy models, the total number of impaired firms 
at time t, 9 tot (t), gives just the evaluation of the total losses due to the propagation 
at this time of the effects of the initial destruction of the single firm in the origin, 
hence it is also a measure of damage spreading: 

N 

e tot {t) = Y,Ht)- (8) 
i=i 

Accordingly, the density can be straightaway obtained as: 

pit) = i - -\e tot {t). (9) 

When taking random delays and/or a random network structure, we study the average 
(Qtot(t)), defined analogously to the average density in Eq. (Q. 
We summarize in [Tab. 2] the different cases that we consider. 



3 Circulant matrix 
3.1 Topology 

Here we look at the network topology corresponding to a braid chain, which is ob- 
tained from a connectivity matrix A (of size N ■ N) with circulant structure. When 
taking periodic boundary conditions for the indexes (i ± N = i), the elements can be 
written as: 

A .. f 1 for j = z — 1, z — 2, . . . , z — n 
lJ \ otherwise, 
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for some integer n < N . This models a directed braid chain, with the same in/out- 
degree n for all the nodes. Each firm is linked to 2n other firms, and needs as 
inputs part of the goods manufactured by the previous n ones (its suppliers), whereas 
the produced outputs are used from the next n ones (its customers). The resulting 
deterministic topology is strongly connected: starting from any node one finds at 
least one directed path along which the signal can propagate to any other node, and 
in fact there are multiple concurrent paths as soon as n > 1. Most of these paths have 
the same time length in the purely deterministic case of equal delays, whereas they 
have usually different time lengths for random delays. In other words, the random 
delays allow to model the asynchrony in the concurrent production paths. We present 
in [Fig. 1] an example for n = 2 in/out- degree. 




Figure 1: The braid chain structure given by a circulant matrix with periodic bound- 
ary conditions in the case of n = 2 in/out-degree. The size of the network is N = 14, 
and the nodes are arbitrarily ordered clockwise. In the deterministic case of equal 
delays also the choice of the origin (the i — 1 node) is completely arbitrary. Note that 
instead, in the presence of stochasticity, the different concurrent paths connecting two 
given nodes do not usually have the same duration in time, although they might have 
the same spatial length in terms of number of steps along the chain. 

When A is a circulant matrix, the equation for the variable of the system ([3]) 
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can be simplified into: 



n 



X 



i(t) = Hi(t) \xi-j{t 




(11) 



b=i 



n 



X 



i(t) = Hi(t) \xi{t 



- Ti.i 



) v Xi-j (t - n 



(12) 



b=i 



for the free models defined by Eq. ([I]) and the forced models defined by Eq. ([2]), 
respectively. The considered systems give an extremely simplified description of the 
input-output structure of real economic network; nevertheless, braid chain structures 
turn out to be an interesting approximation when interpreting each node i, with its 
production corresponding to a whole sector of the economy. 

3.2 The free model with n = 1 and equal delays 

The free model on a braid chain network, with a single input-output connection for 
node, n — 1, and deterministically taken delays all equal to the same time unit To, is 
an example of conservative system of BDEs [211 E2J [231 [31] . Its dynamics is periodic 
from the beginning for all the initial states, without any transient. The model is 
described by the set of equations: 



where we fix initial values Xi(t) = 1 Vz for t G [— To,0), and we represent the initial 
destruction of the single firm on node i — 1, for a duration r c , by the choice of the 
external coefficient: yUi(t) = for t G [0,r c ). The solution displays a wave of nodes 
taking the value zero the one after the other, which propagates periodically across 
the chain, as in the example given in [Fig. 2a]. 

This simple system is analogous to the model introduced in [31] as a first step 
towards a BDE equivalent of hyperbolic partial differential equations. In that work, 
we were discretizing the unidimensional space lattice, and in the analogy the node 
index i plays the role of the discretized coordinate. We studied the evolution starting 
from an initial state where only the variable associated to the lattice point in the 
origin was equal to one, and we found correspondingly a propagating wave of points 
taking the value one in the spatio-temporal pattern. 

In fact, when the duration of the initial perturbation, r c , is smaller than the 
unit of time, To, the role of the external coefficients can be absorbed into the initial 
conditions. The evolution shows in any case no transient and the period of the solution 
is it = Nr . For values of r c multiple of To, the density of firms which maintain their 
production takes the constant value p(t) = 1 — t c /(Nt ), whereas more generally it 
has period r . 

This dynamics implies that, on the one hand, the damage does not spread, but, on 
the other hand, the activity never recovers completely: each impaired firm gets back to 
the unaffected state after the time r c , but contemporaneously its production shortage 
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N. 



(13) 
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Figure 2: Time evolution of the state of the 10 nodes of the free model on a braid 
chain, with n — 1, after an initial perturbation of node i — 1. Each Boolean variable 
is plotted as function of i and it is a piecewise constant function varying between 
and 1. We moreover shift the solutions Xi(t) along the y axis, as a function of the 
node index i, for visualizing the propagating wave in the network. On the left, the 
delays are deterministically taken all equal to To = 1 day, and the duration of the 
initial perturbation is r c = tq/2 = 0.5 day. On the right, we look at a particular 
configuration of random delays, {t^-i}, uniformly distributed between r m j n = 1 day 
and T max =10 days, yielding a period of 47 days. Here the duration of the initial 
perturbation is r c = r m j n = 1 day. 



reaches its customer, with a constant delay To- This unrealistic result is linked to one 
of our simplifying assumptions, namely the discretization of firm production capacity, 
with no possibility for overproduction or production rescheduling. 

3.3 The free model with n = 1 and random delays 

In the free model on a braid chain network with n = 1, the effect of stochasticity in 
the delays can be straightaway worked out explicitly: 

Xiit) = Xi-x(t - r M _i) = 

= Xjit-) y Ti- k ,i-k-l j • (14) 

Correspondingly, one finds that the solution is still periodic from the beginning for 
all the initial states, and the system is therefore once again conservative. 

In the considered case of an initial perturbation of node % — 1 of duration r c , we 
get in particular a periodically propagating wave of nodes taking the value zero the 
one after the other in the spatio-temporal pattern, as in the deterministic case. The 
difference is that the propagation time of the perturbation, from a supplier i — 1 to 
its customer i, is now given by the quenched random variable r^-i, and the period ir 
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of the solution, for a given disorder configuration, is equal to the sum of these delays 
along the whole chain (see [Fig. 2b]), with in average (n) = Nfcj-i) = N(r max + l)/2. 




5 10 15 20 25 30 100 200 300 400 500 

t t 

Figure 3: Time plot of the average fraction of fully active firms, (p{t)), in the free 
model on a braid chain, with n — 1, and stochasticity in the delays. Here r max = 10 
days, as in [Fig. 2b], and we compare the behaviors for network sizes iV = 10, 20, 50 
and 100 with the expected ones (see Appendix A). On the left we present the details 
at short times, and on the right the whole considered time window. The data are 
averaged over M s different random configurations of the delays, a number which is 
taken enough large to give errors of the order of the point size in the plot. 



The density averaged over the disorder, (p(t)), can be computed analytically by 
applying the central limit theorem (see Appendix A). In [Fig. 3], we compare the 
expected behavior, given by Eq. ( )57|) . with the numerical results, for different network 
sizes N: at short times, one can observe the jumps at t = r max ; at long times, there 
are corrections to the expected behavior, since in the considered sums of random 
delays the same variables appear more than once; nevertheless, the mean asymptotic 
value is in agreement with Eq. ( 159]) . 

The case r c ^ r min can be understood within the same framework, and one finds 
that, for large N values, the average density does only depend on the ratio of the 
duration of the initial damage, r c , to the size of the network N itself: 

(p aB (*)>~l-^(fl«(r c = l,t)>. (15) 

Here we are looking at (p a v{t)), assumed to be averaged also over a time window of 
length T min = 1 day: in fact, for r c not multiple of r min , the density (pit)) approx- 
imately oscillates between the expected solution for [r c /r min ] and for [r c /r min ] + 1 
(where [y] means the integer part of y), with period r m j n : this result becomes correct 
in the asymptotic large time limit. 

3.4 The free model with n > 1 and equal delays 

Generally, the BDE systems which describe the free models on a braid chain, with n 
input-output connections, can be reduced to closed sets of n equations; these equa- 
tions involve the products of only n variables, corresponding to nodes in consecutive 
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positions along the chain, but each variable appears more than once, at different 
times. 

In particular, in the deterministic case of delays all equal to To, and taking N 
multiple of n for simplicity, one has: 



Xi(t) 

x i+ i(t) 



Xiit - kr ) 
. ■ x i+ i(t - kr ) 



N-n+1 

x i+n _i(t - kr ) ■ x i+n ^ 2 {t ~ kr ) ■ 

k=N/n 
N-n+2 

x i+n -x(t - kr ) ■ x i+n - 2 (t - kr ) 

k=N/n 
N-n+2 

Xi(t - kr ) 

k=N/n+l 



N-n+2 

x i+n _i(t - kr ) 

k=N/n 
N-n+2 

x i+n _ 2 (t - kr ) ■ ... ■ x i+ x(t - kr ) ■ Xi(t - kr ) 

k=N/n+l 



(16) 



In this expression they are contained all the possible paths along which the infor- 
mation, here the damage due to the initial perturbation, can propagate across the 
network. The reduction is analogous to the one of systems of differential equations 
to systems of higher order equations depending upon a smaller number of variables. 

As soon as the number of connections is n > 1, the dynamics is dissipative [21~ | l22~ | 
|2"3"| 131 j: the asymptotically stable solution is the state in which the whole economy is 
attained by the consequences of the initial damage, Xj = OVz, and it is finally reached 
for any duration of this starting perturbation r c > r (see [Fig. 4]). This result 
can be qualitatively explained by noticing that the present economic models lack of 
flexibility, since the n different inputs to a given firm are considered to be linked by 
AND logical operators. Besides, the system is isolated and the deterministic topology 
of the braid chain is strongly connected. 

We analyze the dynamics of the model starting from a situation in which the 
impaired firms occupy nodes in consecutive positions along the braid chain: let us 
say that in the interval [t,t + 1) (in unity of To) there are 9t t(t) = k impaired firms, 
i.e., Xiit) = for i G [i m (t), ?m(0]> where iuif) = i m (t) + k — 1 and k N. At 
this point, it is helpful to imagine a clock, with the hour-hand marking the position 
of the impaired firm closest to the origin, i m (t), and the minute-hand marking the 
position of the one farthest from the origin, zju(^), at time t. One finds that both the 
hour-hand and the minute-hand move with constant velocities, given respectively by 
1 and n, and therefore (as soon as n > 1), the width of the set of firms which are 
unable to fully produce increases with constant velocity, given by n — 1, or in other 
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Figure 4: Time evolution of the state of the 10 nodes of the free model on a braid 
chain, with n = 2, after an initial perturbation of node 1. We plot Xi(t) as a function 
of t, the solutions corresponding to different node index i being shifted along the y 
axis in order to visualize the spatio-temporal pattern. Here all the delays are equal 
to To = 1 day. On the left, the duration of the initial perturbation is r c = tq = 1 day, 
whereas on the right we take r c = r /2 = 0.5 days. See the text for details. 



words the damage spreads to a total number of n — 1 more firms at each subsequent 
time step. 

In fact, at each time step if — t + 1, the firms in the next n positions after 
will be impaired, since at least one of the stocks that they need is finished, hence 
^m(£ + 1) = Hi{t) + n. The first firm in the sequence is instead the only one which 
recovers its unaffected state, all its suppliers being previously fully active, hence 
i m (t + 1) = i m (t) + 1. Correspondingly, one gets just 9 to t(t') = tot (t) + (n — 1). The 
same argument can be applied again, starting from the new sequence, and so on. 

Following this analysis, the evolution will stop at the time ttrans > when the minute- 
hand reaches the hour-hand with one turn of advantage: the system has attained the 
asymptotically stable steady state a?j = Vi, and the activity can not fully recover 
at all. Notice that the result is unchanged if, in the last time step of the transient, 
the minute-hand surpasses the hour-hand, since Xi = describes a shortage in the 
production of firm i aside from the number of stocks of goods which are not available. 
Moreover, during the transient, each firm can recover at most once. 

Summarizing, 6 to t(t) is linearly increasing with t up to reach the system size N; 
correspondingly, for an initial perturbation destroying the production of one firm for 
a duration r c = To = 1 day, p{t) = 1 — 1/N — (n — l)t/N for t < t tra ns and p(t) = 
for t > t trans , where the length of the transient is straightaway given by: 

N-l 

ttrans — ~T - (17) 

n — 1 

When generalizing to durations of the initial event r c ^ r , it is important to 
distinguish between: 

• t c > t : Here the clock-argument turns out to be valid after the first ~ t c /t 
time steps; therefore, at large times, the damage spreading velocity is constant 
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and equal to n — 1, and the asymptotic stable state is the one where all the 
firms are impaired. 



• t c < Tq\ As shown by data in [Fig. 4b], here one can say that it is the behavior 
of the initially damaged firm, 



that spreads across the network, with constant velocity n — 1; the asymptotic 
solution is periodic of period equal to the unit of time, To, with all the firms 
synchronously impaired only in the first r c part of each period; the density 
p av {t), averaged over the period, is linearly decreasing with time, in this case 
up to the asymptotic constant value t c /tq, reached after a transient of length 
given again by Eq. ( TT71) . 

These results are confirmed by the analysis that we present in Appendix B, where 
9 to t(t) is computed explicitly. 

3.5 The free model with n > 1 and random delays 

We now consider the free model on a braid chain, with multiple n > 1 input-output 
connections, in the presence of stochasticity in the delays. Here the {t^ } are quenched 
independent random variables, taken uniformly distributed in the interval [r m ; n , T max ] 
accordingly to Eq. (T4]), and the unit of time is given by r min , that we choose again 
equal to 1 day. 

As usual, we look in particular at the behavior of the system after the destruction 
of the single firm in the origin for a duration r c . For r c > T max , we expect the damage 
to spread all over the network. In fact, the damage can not spread more slowly than 
in the peculiar case in which all the delays are equal to r max : it follows that the state 
Xi = Vz is asymptotically stable also in the presence of stochasticity, and that it is 
surely reached, after a transient t trans . 

Nevertheless, now the concurrent paths, with the same space length along the 
chain, usually do no longer have the same time length, since in the sums of the involved 
delays appear different random variables. The presence of asynchrony has multiple 
consequences (see [Fig. 5]): the reduced set of equations (Tl6|) does generally contain 
a number of terms increasing with the path space lengths; during the dynamics, there 
is often no more a sequence of consecutive positions along the chain all occupied by 
impaired firms; as soon as r c < r max , each firm can be impaired and recover more 
than once; most intriguingly, because of this reason, the dynamics can be ruled by a 
faster time scale than the average delay. 

In fact, since the delays are independent identically distributed random variables 
along the whole chain, in the limit 1 <C t <C t tra ns, one can argue that the dam- 
age spreads with constant average velocity a, where a m ,i n < a < a max : the lower 
(respectively upper) limit is easily obtained by taking all the delays equal to T min 




for t £ [0, r c ) 
for t G [t c , 1), 



(18) 
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Figure 5: Time evolution of the state of the 10 nodes of the free model on a braid 
chain, with n — 2, after an initial perturbation of node 1. As in the previous cases, we 
plot Xi(t) as a function of t and we shift the solutions corresponding to different node 
index i along the y axis in order to visualize the spatio-temporal pattern. We look 
at a particular configuration of random delays, {t^}, uniformly distributed between 
T m in= 1 day and T max =10 days. The duration of the initial perturbation is r c = r m j n 
on the left, and r c = r max on the right, which yield transients of length 35 and 33, 
respectively. 



(respectively to T max ), hence: 

n — 1 

<a<n-l. (19) 

Tmax 

Nevertheless, as we are going to discuss, the upper limit is usually a definitely better 
approximation to the average velocity than the lower one. 

In order to understand this point, we use once again the analogy with a clock, 
where here the hour-hand marks the position i m (t) of the impaired firm nearest to 
the origin and the minute-hand the position ij^(t) of the impaired firm farthest from 
the origin: the key ingredients are that the network is braid chain structured and 
that the firm is affected by the consequences of the initial event, being unable to fully 
produce, as soon as a single one of the n stocks that it needs is unavailable. Therefore, 
in the long time limit, for 1 <C T max <C n <C N, the average velocity of the hour-hand 
is negligible, (i m (t)) ~ const, and the most of the region between the origin and the 
minute-hand position, is occupied by impaired firms, (9 to t(t)) ~ (hi (t)) —const. 

In other words, the long term dynamics is dominated by the hare, moving more quickly 
than the average runner. 

In this limit, one can moreover evaluate how far does the hare go in one time step, 
of length r m i n = 1, hence the approximated average velocity a* of the signal: from a 
given supplier j, the damage can spread up to the farthest customer i with nj = 1 
delay value; at most, it can move up to i = j + n, with probability V(rj +n j = 1) = 
^/ r max- Generally, the probability that it moves of a distance at most k is obtained 
by taking t^j > 1 for h = j + n, j + n — 1, . . . , j + k + 1 and rj + kj = 1. To get the 
average distance, one has to sum up each possible value k — 0, l...n, multiplied for 
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its corresponding probability, implying that % — j + k is the farthest reachable point; 
hence it follows: 

n u / i \ n ~ k 

«* = ]T 1 -n-(w-l), (20) 

where we used the expansions: 

oo 

fc=0 
oo 

fc=0 

with e = (1 - l/r max ). 

This result is still an underestimation of the effective average velocity in the con- 
sidered limit, since there are corrections of order e n which make a larger, and we 
are moreover neglecting that the signal can go still faster on more than one time 
step. However, notice that, in a large region of the parameter values, one finds 
a* > {n — l)(r max + l)/r max , thereby confirming that the dynamics is usually domi- 
nated by a faster time scale than the one corresponding to the average delay. 

A different approach is worked out in Appendix B, where the average number of 
impaired firms is explicitly computed as a function of the probability for the signal to 
have propagated of k position in t time steps: this analysis confirms the expectation 
of a linear behavior on a large time window, and it allows to predict the effective 
slope. 

We present in [Fig. 6] our results on the average density of fully active firms, 
(p{t)), which is in fact linearly decreasing with time on a large window. The obtained 
behavior is well in agreement with our expectation, since the data lie between the 
lower limit p m in{t) ~ 1 — {n — l)t/N, corresponding to the peculiar case in which 
all the delays are equal to T min , and the upper limit p max (t) ~ 1 — a*t/N, evaluated 
using the hare argument. Moreover, the approach worked out in Appendix B allows 
to get results on (p(t)) which are nearly indistinguishable from numerical data. 

Here we also consider delays depending only upon the customers, = r(i) Wj, 
or only upon the suppliers, T^j = r(j) Vz, where {r(i)} (respectively {r(j)}) are, 
as usual, independent identically distributed random variables, which follow the law 
Eq. (jlj): we get the same, slightly slower, average signal velocity in these last two 
cases. This result can be qualitatively explained by the observation that here there 
is a smaller number of propagation paths of different time lengths. 

We found the same kind of almost linear decay, with velocity quite correctly ap- 
proximated by a*, also for different values of the model parameters. We moreover 
checked that, as soon as r c > r min , the average density shows a dependence on r c 
only on the first ~ r c time steps. Finally, we studied the behavior of the densities 
p(t) obtained from single different disorder configurations: this quantity does usually 
slightly fluctuate around the average value. In particular, the probability distribu- 
tions of the transient length tt ra ns is Gaussian-shaped, becoming rapidly peaked for 
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Figure 6: The average density of fully active firms, (p(t)), as a function of time, in 
the free model on a braid chain, after an initial perturbation of node 1 of duration 
i~ c — T-min- Here the network size is N = 10000, and the number of input-output 
connections is n = 20; the average is taken over J\f s = 100 different configurations of 
the disorder. We study delays depending only upon the customers, Tij = r(i) Vj, 
only upon the suppliers, Tij = r(j) Vz, and upon both of them, nj = r(i,j); in all 
of the cases they are independent random variables, identically uniformly distributed 
according to Eq. (TJJ, with T min = 1 day and r max = 10 days. We compare the data 
with the lower (a = n — 1), with the upper (a = a* = n — {T max — 1)) limits on the 
behavior, and with the prediction obtained by the approach explained in Appendix B. 
See the text for details. 
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increasing network size N: this means that most of the configurations of the delays 
gives a solution which reaches the asymptotic one, Xj = OVi, in about the same time, 
of order t tra ns ~ a*/N. In other words, the density and t tra ns are self-averaging quan- 
tities, whose values are independent on the particular disorder configuration in the 
N — > oo limit. 

Anyway, from the point of view of the dependence of the asymptotic solution 
upon the disorder, it is interesting to notice that the steady state x« = Wi is surely 
reached only for r c > T max : otherwise, there are configurations of the delays in which 
the smallest variable is minij{Tjj} = r* > r c and, for the same argument that we 
discuss in the deterministic case when r c < tq, the asymptotic solutions in these cases 
turn out to be periodic of period r*, with all the firms simultaneously impaired in 
the first part of length r c of the period. 

Nevertheless, since the probability of obtaining a configuration of the delays in 
which minjj jrjj} = r*, V T *, approaches rapidly zero for increasing N values, in the 
limit of large network size N, periodic asymptotic solutions are usually found only 
when the duration of the initial damage r c is smaller than the smallest possible nearest 
neighbor time length, r m j n = 1. In fact, this probability V T * can be straightaway 
evaluated: 

/ t*-1\ n2 

V T * = 1 for r min <r*< r max , r* E K (22) 



3.6 The forced models with n = 1 

We now turn onto the study of the forced models on a braid chain, defined accordingly 
to Eq. (|12p . When the in/out degree is n — 1, and the delays are deterministically 
chosen all equal to r , the BDE system Q can be reduced to the equation: 

= { Yl Wi ~ k [ t ~( k + i)^] } V Xi {t - Nt q ), (23) 

where also the sum in the parenthesis is to be interpreted in the sense of the OR 
Boolean operator, V. For randomly distributed delays the equation is slightly more 
complex: 



Xi{t) 



'JV-l / k 

.k=0 \ j=0 




(24) 



since the time lengths of the paths are here given by the sums of the corresponding 
random variables. 

One can distinguish between: 

• A duration of the initial damage r c smaller than the smallest nearest neighbor 
propagation path, which means: 

— Tc < r o i n the deterministic model; 
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Figure 7: Time evolution of the state of the 10 nodes of forced models on a braid 
chain, with n — 1, after an initial perturbation of node 1. We plot Xi(t) as a function 
of t and we shift the solutions corresponding to different node index i along the y 
axis in order to visualize the spatio-temporal pattern. On the left we present data for 
the deterministic case of delays all equal to To = 1 day, whereas on the right we look 
at a particular configuration of random delays {t^}, uniformly distributed between 
T m in = 1 day and r max = 10 days, with min, r^j-i = t* = r m %n- The duration of the 
initial damage is taken to be t c = 2tq = 2 days and t c = (rj )3 -) = 5.5 days, respectively. 
See the text for details. 



— t c < mini j {r» j-} = t* in the stochastic one, where the probability to have 
t* > T. m in is given by Eq. fl22} , and becomes negligible in the limit of large 
system size. 

Here the systems are conservative and the solutions are the same as in the 
corresponding free models (see [Fig. 2]). In particular, they are periodic from 
the beginning, with a wave of impaired firms which propagates in the spatio- 
temporal pattern. The periods are given by the sums of the delays along the 
whole chain, hence: 

— ii = Ntq in the deterministic case; 

— 7r = J2f=i r i,i-i) with (tt) = N(T max + l)/2, in the stochastic one. 

• A duration of the initial damage t c larger than the smallest nearest neighbor 
propagation path (t c > To or t c > t* in the deterministic or stochastic model, 
respectively). As shown in [Fig. 7], here the solutions are not periodic from the 
very beginning. In detail, when the firm in position i is reached by the wave of 
damage, its activity is affected for a duration: 

— no longer than To in the deterministic case, 

— no longer than r^-i in the stochastic one. 

The asymptotic solutions are surely reached after a first cycle in which each 
firm is impaired and recovers once. Moreover: 
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Figure 8: On the left, we present the time evolution of the state of the 10 nodes of the 
deterministic forced model on a braid chain, with n — 2, after an initial perturbation 
of node 1. As in the previous considered cases, we plot Xi(t) as a function of t by 
shifting the solutions corresponding to different node index i along the y axis. On the 
right, we show the density of fully active firms as a function of time in the same model, 
for a network size N = 10000, and different values of the input-output connections 
n. In both of the cases we take a duration of the initial damage r c = r = 1 day 

— they are the same as for r c = tq in the deterministic model, 

— they are the same as for r c = r* (usually equal to r m j n in the large size 
limit), in the stochastic one. 

3.7 The deterministic forced model with n > 1 

The behavior of the forced models on a braid chain, with in/out-degree n > 1, can 
be understood in detail when the delays are deterministically taken all equal to To- 
Here, Eq. ffl2|) simplifies into: 



Therefore, the firm i maintains its production, as usual, if all its suppliers were fully 
active at the previous time step, but also if some of them were impaired and i itself 
was already impaired: the production is recovered after one time step of length r . 

Let us take for simplicity N multiple of n, and let us start by considering a 
duration of the initial damage r c > tq: then, after the first time step, there are n firms 
in consecutive positions along the chain which activity is simultaneously impaired; at 
the subsequent time step, these firms recover but the damage propagates to the next 
n ones. 

In other words, in this model, as soon as t > 1, both the hour-hand and the 
minute-hand move with the same constant velocity n: the damage does not spread, 
but the activity never recovers completely. Correspondingly, the asymptotic solution 




(25) 
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is periodic of period N/n and the density is just p(t) = 1 — n/N, as shown by the 
data presented in [Fig. 8]. 

Though here there is no complete breakdown of the economy, this simple case 
makes evident that, in the considered models, the presence of more connections does 
lead to a less favorable final outcome, with more impaired firms. This result can 
be explained by the fact that more connections do not lead to risk sharing, since 
having one impaired supplier is enough to stop a firm production. This assumption 
amounts to say that each supplier of a firm provides a different type of goods or 
services, and that one supplier cannot compensate for the loss of another supplier. In 
such a situation, a firm that is connected to more suppliers has higher risks of being 
indirectly affected by a shock. 

When r c < r , one still finds a propagating wave of n impaired firms in consecutive 
position along the chain; nevertheless, here their activity is simultaneously impaired 
only in the first part of length r c of the time step. The situation is similar to the one 
encountered in the same case in the free models with n — 1, i.e. it is the behavior of 
the initially attained firm (see Eq. f fT8|) ) which propagates along the network. 

3.8 The stochastic forced model with n > 1 
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Figure 9: Time evolution of the state of the 10 nodes of the stochastic forced model 
on a braid chain, with n = 2, after an initial perturbation of node 1. We plot Xj(t) as 
a function of t by shifting the solutions corresponding to different node index i along 
the y axis in order to visualize the spatio-temporal pattern. As usual, we look at a 
particular configuration of random delays, {t^}, uniformly distributed between r m j n = 
1 day and T max =10 days. The duration of the initial perturbation is r c = r m j n . Notice 
that the solution does not display any periodicity in the considered time window. 



In the presence of stochasticity in the delays, the dynamics of the forced models 
on a braid chain, with n > 1 input-output connections, turns out instead to be quite 
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complex. Here we take as usual random independent delays, uniformly identically 
distributed according to Eq. (jlj), with r max = 10, in the unit of time given by r min = 1 
day. The solution for a typical disorder configuration, after an initial damage of 
duration r c = T m j n , is shown in [Fig. 9]: intriguingly, despite of the relatively small 
network size considered, N = 10, and of the small n = 2 value, one does not observe 
any periodicity in the considered time window. 

Since the delays are integer multiple of T min , the asymptotic solution of the system 
([3]), for a given disorder configuration, is surely constant or periodic because of general 
results on BDEs [22], [23], [3T] . Nevertheless, in the same works it is shown that, for 
irrationally related delays, BDEs can have solutions of increasing complexity, which 
display a number of jumps in time intervals of the same length increasing with time; 
this peculiar behavior is observed in particular in the case of linear systems such as 
the simple example: 

x{t)=x{t-r)\/x{t-l). (26) 

Here, r G (0, 1) is irrational, and v means the XOR Boolean operator, hence x(t) is 
true (x(t) = 1) only if a single one among the two variables x(t — r) and x(t — 1) is 
true. It can moreover be shown that these aperiodic solutions can be approximated 
with the desired accuracy, for increasingly long times, by the periodic solutions of 
nearby BDEs systems with rationally related delays, which approximate better and 
better the irrational ones. 

Though a careful analysis would be necessary in order to extend these results to 
the present case, we notice that Eq. f[T2l) . describing the model, is equivalent to: 

Xi(t) = fii(t) i^y~)xi(tjj) V Xi-j{tij) V fiiiUj) ■ j , (27) 

where we label ty = t — r^j-j. The equivalence is due to the Boolean relation a V b = 
a Ab = a y b y (a • b) , and it makes clear that the considered system is partially 
linear. This suggests that the long periods observed in the behavior of the solutions 
can be related to the known results on the existence of small, partially linear BDEs, 
with rationally related delays, which approximate the increasingly complex solutions 
of systems with irrationally related delays. 

In order to better characterize our numerical findings, we find interesting to study 
various quantities: 

• The length of the transient, T trans , rigorously defined just as the time elapsed 
before that periodicity settles in. In the previously considered stochastic free 
models, T trans is the same as the time t trans , which the average density takes 
to reach the asymptotic zero value. Instead, as we are going to discuss, in the 
present case one finds T trans 3> t trans , and it is therefore important to distinguish 
between the different meanings of the two transients. 

• The period n of the asymptotic solution itself, possibly equal to zero if it is 
constant. 
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The period-averaged asymptotic density of fully active firms, p^, defined as 



1 /"Ttrans+T 

= — j p{t)dt. 
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Figure 10: Data on the forced model on a braid chain with random delays, uniformly 
distributed between r m j n = 1 day and T max = 10 days. Here we study the solutions of 
small systems of N ~ 10 equations with n — 2, after an initial perturbation of node 
1 of duration r c = r min . The results are obtained from Af s = 1000 different disorder 
configurations. On the left we plot the probability distribution of the transient, 
P(T trans ), in the center the probability distribution of the period, -P(vr), and on the 
right the probability distribution of the period-averaged asymptotic density of fully 
active firms, P(p OQ ). Both P{T trans ) and P(tt) are presented in log-log scale, in order 
to make evident that they are not negligible also at very large values. See the text 
for details. 



These quantities can be straightaway exactly measured up to relatively small sys- 
tem sizes N ~ 12. We present in [Fig. 10] our results on their probability distribu- 
tions, P(T trans ), P(tt), and P(poo) respectively, as obtained by considering Af s = 1000 
different configurations of the random delays. 

The most striking feature displayed by the data is that both the transient, T tra ns, 
and the period, ir, of the solution can increase very fast with the network size N. In 
order to make it evident, we have plotted the resulting probability distributions in 
log-log scale (see [Fig. 10a] and [Fig. 10b]): for instance, P(T tra ns) is not negligible 
for values of the transient as large as T tra ns ~ 10 6 (in units of T m i n ), in systems of 
N = 12 nodes. In detail, we found that the average transient diverges exponentially 
with N, (Ttrans) oc exp(constN) , and that, at least for the considered N values, the 
average period (it) displays a similar behavior. 

The data on the asymptotic density of fully active firms, poo, averaged over the 
disorder dependent period n, are shown in [Fig. 10c]: its probability distribution, 
P(poo), turns out to be bell-shaped, and it becomes more peaked around the mean 
value for increasing network size N, therefore suggesting that the quite complex 
dynamics does not imply the unpredictability of the behavior of macroscopic intensive 
quantities, at least in the large network size limit. 

In fact, whereas to predict the details of the solutions for given disorder configu- 
rations of the random delays seems to be a quite hard task, one can use an approach 
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similar to the one of Section 3.5 to obtain the expected behavior of various average 
quantities. 

For simplicity, we limit the analysis to the case in which the duration of the 
initial damage is not smaller than the smallest nearest neighbor propagation path, 
t~c > T m i n . In the present extension of the clock argument, we assume that the 
minute-hand marks the position iiaif) of the firm farthest from the origin, i.e. from 
the node where the initial damage acted, which has already been impaired at least for 
a duration T min at time t. We recall that the network is braid chain structured and 
that the production of a given firm, apart from the initially damaged one, is impaired 
for the first time as soon as one of the n stocks that it needs is unavailable. For these 
reasons, in the limit 1 <C r max <C n <C N, at large t, the signal is propagating across 
the chain with velocity still well approximated by Eq. (120]) : 

(i M (t)) ~ aft =[n- (r max - l)]i. (29) 

The main difference with the previous stochastic free model is that here, because 
of the external rescue inputs, only a fraction s* of the firms between the origin and 
(*m(£)) remains in average impaired during the time step of length r min . Since the 
delays are taken independently and identically distributed along the whole chain, at 
large enough times this fraction is constant, and the damage spreads, up to invade 
the whole network, according to: 

(M*)> ~ s * a *t = s*[n- (w - l)]t, for 1 < t<t trans . (30) 

Correspondingly, the average density of fully active firms is again linearly decreasing 
with t, up to the time t trans , at which it reaches a nearly constant asymptotic value: 

«t)>~(;- ( f /jv)a * f ' t j** 5 *— pi) 

[ t — S , IOr t > ttrans- 

Hence, the effective transient, t tra ns, relates to dynamics of (p(t)); in fact, it is the 
time at which the minute-hand reaches again the origin after a whole tour, and it is 
therefore of the same order as the time before that the density became zero, since the 
system attained the asymptotically stable configuration, Xj = Wi, in the stochastic 
free model: 

N N 

ttrans ~ — — 7 7T- (32) 

a* U- (Tmax ~ 1) 

Therefore, whereas in the free model one has t tra ns = Tt rans , in the forced case it is 
important to distinguish between the transient T tra ns, rigorously defined as the time 
elapsed before that periodicity settles in, that, as we showed numerically, is in average 
exponentially increasing with the network size N, and the definitely smaller effective 
transient t t ran s , which can be more generally defined as the time at which macroscopic 
observables attain nearly constant average values. 

We present in [Fig. 11] the observed behavior of the density of fully active firm 
as a function of time, after an initial perturbation of node 1 of duration r c = r m j n , 
for a large network size N = 10000 and in/out-degree n = 20. In agreement with 
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Figure 11: The density of fully active firms p(t) as a function of time, in the forced 
model on a braid chain with random delays, uniformly distributed between r min = 1 
day and T max = 10 days, after an initial perturbation of node 1 of duration r c = r m j n . 
The network size is N = 10000 and the in/out-degree is n = 20. The results are for 
a single (typical) disorder configuration. We study delays depending only upon the 
customers (r^ = r(i) Vj), delays depending only upon the suppliers (r^ = r(j) Vi), 
and delays depending upon both of them (ry = r(i,j)). See the text for details. 



the picture emerging from the previous discussion, we find a nearly linear decay on a 
quite large time window, followed by small oscillations around the asymptotic value, 
which is reached in a time t trans of the same order of the estimation of the effective 
transient, t tr ans ~ 900 (in units of r m j n ), that one gets from Eq. fl32|) . From this 
point of view, we also recall that our evaluation of a* is an underestimation of the 
damage spreading velocity, hence the corresponding value of t trans is to be considered 
an upper limit. 

The data are for a single, typical, disorder configuration, in order to make evident 
that the fluctuations of p(t) around the average are small, as it can be expected for 
a macroscopic intensive quantity: we checked in particular that both the effective 
transient length, ttrans, and the asymptotic mean value of the density, 1 — s*, are 
usually the same for different choices of the random variables. 

Besides considering delays depending both upon the customers and the suppliers, 
we present, in the same [Fig. 11], data on the system with delays depending only 
upon the customers, nj = r(i) Vj, or only upon the suppliers, r^j = r(j) Vz: the 
last two cases are characterized by compatible values of the effective transient, ttrans, 
which is slightly larger than in the first one. This can be explaining by noticing 
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that, as in the stochastic free model, when Tjj depend on both of the indexes, there 
are more concurrent paths of different time lengths in the system, hence the average 
signal propagation velocity is higher. 

The nearly constant asymptotic mean value of the density, 1 — s*, is instead 
definitely larger when the delays depend only upon the customers, whereas the other 
two considered cases give very close results. In fact, if r^j = r{i) Vj, Eq. ffT2]) becomes: 

Xi (t) = Vi(t)xi[t - r(i)] V jlJzi-^ - r(i)] J . (33) 

Since the activity of firm i is impaired as soon as a single one of the stocks of products 
that it needs is lacking, the fact that the delays in production are assumed to be 
independent on the particular customer implies that there are much less combinations 
of the delays which can result in shortening the production of a given firm, i.e. the 
average fraction s* of contemporaneously impaired firms, in the region reached by the 
spreading of the initial perturbation, is smaller. 

Moreover, in this particular case, one can make a further step in the analysis: in 
fact, Eq. (133]) implies that, after the transient T tra ns, when averaging the dynamics on 
an enough large time window, each firm is roughly impaired for one half of the time; 
for large N 3> n, this turns out to be equivalent to say that there is in average one 
half of impaired firm at each time step, i.e. s* ~ 0.5. This is just the result on the 
mean asymptotic value of the density observed in [Fig. 11], which we checked to be 
independent upon the considered n. We also found that the solutions for ry = r(i) Wj 
usually display a more regular behavior than in the other situations. 



4 Random graph 
4.1 Topology 

We are interested in the topology induced by a directed random graph, which is 
obtained when the elements of the connectivity matrix A are chosen accordingly to 
Eq. 05]): this is a quite well known random structure [331 El], hence we start by 
summing up some of the main results. 

Random graphs have been initially introduced by Erdos and Renyi about 50 years 
ago [35] [36] EZ] , and they have been extensively studied more recently, when also a 
number of related models have been considered [38] |39] HQ] HI] . One gets an undirected 
Erdos- Renyi random graph, G(N,p), by taking the edges, connecting each possible 
pair of the N nodes, independently identically distributed with probability p. 
In our notation, the matrix A is symmetric, since the event Ay = 1 implies the event 
Aji = 1, and vice- versa. 

The total number of pairs of nodes is N(N — l)/2, and each edge contribute to the 
degree of 2 nodes, therefore the average degree is straightaway z = (k) = (N — l)p ~ 
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Np. More in detail, the probability distribution of the degree k is binomial: 

V{k) =( N k 1 ) ~P) N ^ k - tf-", (34) 

and it converges to a Poisson distribution with average value z in the large N limit. 
One moreover defines the connected components S c of the graph as the ensembles 
of nodes such that, from each node i G S c , there is at least one path, across nodes 
belonging to the component itself, which reach each other possible node j in the same 
connected component: 

e S c ^dj : A ihl ■ A hlfl2 ■ ... ■ A hd ^ hi,h 2 , .-,hi e S c . (35) 

Therefore, the average size of connected components S c can be evaluated by start- 
ing from a randomly chosen node, and then looking at the number of its first neighbors 
Zx, of its second neighbors z 2 , and so on. For a Poisson distribution, which is a quite 
peculiar case, one simply has: 

~2 

hence one obtains: 

S c 

Correspondingly, one finds that, in the N — > oo limit, (S c ) diverges for z /* z c = 
1. Above this "critical" value it appears a giant connected component, S gc , which 
contains a finite fraction of the nodes, S gc = s gc N, such as it is usually observed in 
real networks. This "phase transition" was already enhanced in the pioneering papers 
by Erdos and Renyi [351 ESI [37], and was subsequently studied in great detail both 
from the mathematical point of view [38J , and from the physical point of view [10] . 

A simple argument for evaluating s gc for z > z c runs as follows [531 H2]: us 
assume that they still exist finite size connected components for z > z c ; then, it is 
s gc = l — v(z), where v(z) is the probability for a random node of being in a finite size 
connected component. Therefore, v(z) is the probability that all the nodes which can 
be reached from a random one is finite, but each node generates a Poisson branching 
process (the number of its first neighbors is Poisson distributed), hence v(z) must be 
the solution of the transcendental equation: 

k=0 



oo u 



k=0 

OO y. 



k=0 



Z2 / Z 2 
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Zl \z 1 
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Zi = z 



(36) 
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for z < 1. 
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This result can also be derived in the framework of the probability generating func- 
tions jl3] , which is well suitable for applications to more general degree distributions 
and to directed network, as recalled in Appendix C. 

In fact, the directed random graph that we are considering, D(N,p), is a simple 
generalization of the Erdos-Renyi model: each directed link is chosen independently 
with the same probability p among the N(N — 1) possible ones. Hence the average 
in/out-degree z = {k in ) = (k out ) is again given by z = (N — l)p ~ Np, and V(k in ) = 
'P(k out ) are still described by Eq. flMl) . Notice that, here, the average in/out-degree z 
corresponds to the connectivity n in the deterministic braid chain structure previously 
considered, and furthermore that the total average degree of the node is 2z, i.e., if we 
were to transform the directed random graph in an undirected one, by interpreting 
each link as an edge, we would get an RG with average degree equal to twice the 
average in/out-degree of the starting DRG. 




Figure 12: A possible realization of a random graph (on the left) and of a directed 
random graph (on the right) with the same small network size of N = 12 nodes, and 
the same p ~ 1/N (hence z ~ 1) value. The number of the edges in the RG is one 
half of the number of the directed links in the DRG. 

We present in [Fig. 12] an example of the kind of structures which can be obtained: 
we compare in particular an RG and a DRG of equal small size N = 12, with the 
same p ~ 1/N (hence z ~ 1) value. This figure makes evident the completely different 
topology with respect to the brain chain, where for the same connectivity value n = 1 
one would find a single connected component with directed links i — > i + 1 (see [Fig. 
1]). Instead, in the shown example, looking at the case of the undirected random 
graph for simplicity, we observe two isolated nodes, three connected component of 
size 2, and one connected component of size 4. 

The figure also makes evident that in the directed case [121 H3j HU 05] , since the 
existence of a path connecting i to j does not usually imply the existence of a path 
connecting j to i, for each given node one can define: 

• the out- component, which is the set of all nodes that can be reached from it; 

• the in-component, which is the set of all nodes from which it can be reached; 
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• the strongly connected component, which is the set of all the nodes that can be 
reached from it and from which it can be reached, hence the intersection of the 
in-component and the out- component; 

• the weakly connected component, which is the set of all the nodes that can 
be reached from it or from which it can be reached, hence the union of the 
in-component and the out-component; this corresponds to the connected com- 
ponent of the node in the graph obtained by disregarding the directionality. 




Figure 13: A sketch of the bow tie topology structure characteristic of directed net- 
work [H]: with respect to the notation used in the text, the two bows correspond to 
the giant components X \ S sc and O \ S sc respectively, whereas the tie represents the 
giant strongly connected component S sc = Xfl O. Here we consider the most general 
case in which W = X U O U T, where T contains in particular the paths linking the 
two bows without passing across the tie, which have been observed in real network 
such as the web [15]. 

The analogous of the "phase transition" in an undirected random graph can now 
be characterized by the possible formation: of a giant in-component I, containing 
/ = siN nodes; of a giant out- component, O, containing O = sqN nodes; of a giant 
strongly connected component, S sc = Xfl O, containing S sc = s sc N nodes; of a giant 
weakly connected component W, containing W = s w N nodes. In fact, one usually 
expects to observe two different transitions, i. e. two different abrupt changes in the 
properties of the system in the large iV limit: at the lower average in/out-degree z™, 
corresponding to the critical average degree in the undirected graph, where it appears 
the giant weakly connected component, and at the higher average in/out-degree zf, 
which gives the effective transition in the directed network, where the other giant 
components appear contemporaneously. 
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For z > z£ the resulting picture is captured by the bow tie structure sketched in 
[Fig. 13], which has been observed in many different real networks [4"5ll4"6] : notice that 
the appearance of the giant in-component (respectively out- component) corresponds 
to the divergence of the number of nodes which can be reached from a given one 
(respectively from which a given one can be reached), i.e. a node is in the giant 
in-component if its out-component diverges, and vice-versa. Moreover, above the 
transition, it can be necessary to introduce one more giant component in order to 
fully characterize the topology, since one can have Wd(IU0): this is in particular 
the case when there are directed paths between X and O which do not pass across the 
strongly connected component S sc , as it has been found in the structure of the web 
[35] , which can be seen as a directed graph characterized by power law in/out-degree 
distributions. 

The formalism of the probability generating functions for directed random graphs 
[121 111] turns out to be particularly simple when the in/out-degree distribution fac- 
torizes, Vk in ,k out = ^W^wt) since the in and out components are independent in the 
large network size N limit, so that one has simply s sc = s/s<> F° r the considered case 
of V(ki n ) = V(k out ) given by Eq. (I34p . the transition does still occur at the critical 
average in/out-degree z% — 1, and the giant in-component and out-component have 
the same size: 

si = s = 1 - v(z), (39) 

where v(z) is once more given by the solution of Eq. (138]) . Then, s sc = (1 — v(z)) 2 , 
whereas the fraction of nodes in I \ S sc , and equivalently in O \ S sc , is equal to 
v(z)(l — v(z)): the three regions in the bow tie have roughly the same size in the 
thermodynamic limit, as it has been observed in real networks [15], for v ~ 1/2, 
hence z ~ log 4. Moreover, one simply has z™ = zf/2: since when v{z) is solution of 
Eq. (I38I) one has v(2z) = v 2 (z), it follows that, in the region z > zf, s w = 1 — v 2 (z) = 
Si + sq — s sc , which confirms that here W = X U O. 




Figure 14: The local tree-like structure characteristic of Erdos-Renyi random graphs: 
on the left we consider the damage spreading on a directed random graph D(N,p), 
with z = (k ou t) = (ki n ) = 3, and on the right the one in the undirected network with 
the same p value, G(N,p). 

The other key ingredient of the topology of Erdos-Renyi random graph is that 
their structure, for not too large p-values, is tree-like: this is rigorously correct for 
the finite size connected components, both below and above the critical point, since 
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it can be shown that here the probability of closed loops approaches zero in the limit 
N — > oo, but it is locally valid also within the giant connected components, where 
the closed loops can be neglected in a first approximation. We compare in [Fig. 14] 
these kind of structures for a DRG and for an RG with the same z = 3 value: notice 
that in the first case one has Z2 = (ki n ){k out ) = z 2 , where in the second case one gets 
the same results, z<i = z 2 , from Eq. ( 136]) . In other words, in random directed network 
the probability to have both the link i — > j and the link j — > i is ~ p 2 , and it is 
therefore negligible in the region p< 1 that we are studying, whereas in undirected 
random graphs we have to consider the edges emerging from the node apart the one 
along which the signal arrived, in order to make a meaningful computation. 

As we are going to discuss, this difference is evident in the results on the damage 
spreading at short time in the BDE models on random graphs, that can be predicted 
quite accurately on the basis of the local tree-like topology. Moreover, when consid- 
ering randomly distributed delays, though also here two paths with the same space 
length along the network do not usually have the same time length, because of the 
local tree-like topology, the number of different concurrent paths connecting two given 
nodes does not increase as fast with the average in/out-degree as in the previously 
considered deterministic braid chain. 

Finally, we note that on the behavior of D(N,p) there are detailed rigorous math- 
ematical results [33] (see also [3lj and references therein), which would require a 
definitely deeper analysis for being summarized. 

4.2 The free models with equal delays 

We introduce the free models, defined by Eq. (TjQ), with all the delays equal to the 
same time unit To = 1 day, on the directed random graph, D(N,p), and we study as 
usual the consequences of a damage initially destroying a single firm for a duration 
t c . For simplicity, we limit the analysis to the case r c = r . 

In the presence of stochasticity in the network structure, the BDE system ([3]) can 
no more be reduced, and one has in principle to solve a set of iV equations in N 
variables, which can be rewritten in the form: 

Xi(t) =fa(t) JJ Xj(t-T Q ) i = l,...,N. (40) 

Here the product, which means as before the AND Boolean operators A, runs over all 
the values of the index j labeling firms which are suppliers of i, hence belonging to 
the set of nodes X{i) = {j : ^ 0}, which corresponds to the definition of the node 
in-component T{i); the in-degree of the node is the total number of suppliers of firm i, 
i.e. = k in {i). Analogously, the out-component of the node, 0{i) = {j : ^ 0}, 
corresponds to the firms which are customers of i, which total number is k out (i). 

Therefore, Eq. (T4T)]) makes evident that the properties of the BDE system strongly 
depend on the probability distributions of the in/out- degree, which in the consid- 
ered case are independent, and described by Eq. (13~4"|) . with mean values (k in (i)) = 
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(koutii)) = (N — l)p = z. Hence, for increasing z values, the properties of the solu- 
tions will reflect the appearance of giant in/out-connected components in the graph, 
and one will find different kinds of solution for average in/out-degree lower or higher 
than the critical value z~ — z c — 1. 

In detail, the damage spreading in this simple BDE model can be understood with 
the same argument used for evaluating the average size of the connected components 
in the graph: at t — 1, the signal propagates from the node i, occupied by the initially 
damaged firm, to its first neighbors, which average number is z\ — z; at t — 2, it 
reaches the second neighbors, which average number is z% = z 2 , and so on. From 
Eq. ( 1361) it follows that at time t the damage did spread up to reach in average z l 
nodes, hence: 

(6 tot (t)) ~ z* for t < log N/ log z, (DRG). (41) 

This implies that the average number of firms which production is impaired increases 
with time only if z > z^ = 1, i.e. only if the graph is above its transition point. 

The argument does make use of the local tree-like structure, since we are implicitly 
assuming that the probability for two node reached at a given time step of being 
themselves connected by a different path, i. e. the probability of closed loops, can be 
neglected. In fact, this argument breaks down roughly at the time where the signal 
propagated to the whole connected component to which the initial node i belongs, 
and one can find different situations: 

• i is in a connected component containing a small number of nodes, S c = 0(1); in 
this case one expects that there are no loops, hence the system can completely 
recover from the initial damage, i.e. 6 to t(t) = in the large time limit; a ran- 
domly chosen initial node i will belong almost surely to a connected component 
containing a finite number of nodes, in the N — > oo limit, for z < z d , and with 
probability v(z) for z > z d c \ 

• the graph is above the critical point and i belong to a giant connected compo- 
nent; because of the presence of closed loops, here a finite fraction of the firms 
will be impaired in the asymptotic solution, and the economic network does 
never recover completely from the initial destruction event. 

Such results can be explained by noticing that, in our model, there are firms which 
have no "clients" in the network. This surprising possibility simply arises from the 
fact that some firms have only final customers as clients, not other firms. For these 
"final-demand" firms, being unable to produce does not have any consequences on 
the rest of the productive system, only on household well-being. Since z measures 
the average number of clients which are in the network, in the limit of a large total 
number N of firms, the phase-transition just corresponds to the fact that, for z < 1, 
the most of the firms are in this peculiar situation, or have a few customer firms 
which are themselves in this peculiar situation, and so on, so that the initial damage 
usually does not propagate; conversely, as soon as z is larger than one, a finite fraction 
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1 —v{z) of the firms have clients in the network, which have themselves other clients 
in the network, and so on. 

To make more quantitative the analysis, in the large N limit, one can say that 
in directed random network above the transition point, in order to observe damage 
spreading up to reach a finite fraction of the network, the initially attained firm has 
to belong to the giant in-component, i G I: this occurs with probability 1 — v(z). 
Moreover, in the same limit, one expects that finally all the firms in the giant out- 
component O will be impaired, which means a fraction so = 1 — v(z) of the whole 
network. Correspondingly, we find: 

{9tot{t)) ^{[l-v(z)) 2 N for t>tl,iV»l (42) 

where the time that (O to t(t)) takes before reaching the asymptotic constant average 
value is quite short also for large N values, since it is of order t trans ~ log([l — 
v(z)]N)/ log z. 

Interestingly, [1 — v(z)] 2 is also the fraction of nodes which are in the giant strongly 
connected component, S sc , of the network; nevertheless, one should note the different 
meaning of the two quantities: when the initial node belongs to X, the total number 
of firms which activity is finally impaired is [1 — v(z)\N = So, definitely larger than 
S sc for small v(z) values. This point can be better understood by looking at the 
Erdos Renyi undirected random graphs with the same average degree z > 1: here, 
the initially damaged firm i has to belong to the giant connected component S gc , 
and the firms which activity is finally impaired are the ones in S gc as well, hence we 
find again that {6 tot {£)) approaches [1 — v(z)] 2 N at large times, though it is clearly 
S sc =l- v(z). 

Summarizing, both in the DRG and in the RG, the average asymptotic density of 
fully active firms, when the networks is above the transition point, turns out to be: 

(Pasym(z)) = 1 - [1 - v{z)f = 2v(z) - v\z), (43) 

where we recall that v(z) is the solution of Eq, 038p ). Instead, one important difference 
between the behavior of the considered free BDE systems with equal delays on directed 
and undirected random graphs concerns the short time dynamics: as it is shown in 
[Fig. 14] , in the undirected case the signal does also propagate back to the node from 
which it is arrived and, correspondingly, to correctly describe the behavior of {O to t{t)), 
Eq. (J4TJ) should be replaced by: 

[t/2] 

(M*)) = J2 zt ~ 21 for 1 < log ^/ log z, (RG), (44) 

1=0 

where as usual we denote [y] the largest integer smaller than y. 

We present in [Fig. 15] the results of numerical simulations of the BDE free model 
with equal delays on random graphs with a large number of nodes N = 5000 -v- 10000, 
for different values of z > 1, above the transition point. On the left, we compare the 
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Figure 15: Data on free models on random graphs with a large number of nodes 
N 1 and deterministically chosen delays all equal to To = 1 day. The obtained 
results are averaged over at least J\f s = 200 different configurations of the links. On 
the left, we present the short time behavior of the average total number of impaired 
firms, (9 to t(t)), both in the directed random graph, D(N,p), with average in/out- 
degree z ~ Np = 1.5 and in the undirected Erdos-Renyi one, G(N,p), with the same 
p value; here N = 5000. On the right, we plot (9 to t{t)) in the directed random graph, 
D(N,p), for different z ~ pN values above the transition point and N = 10000. The 
obtained curves are compared with the corresponding expected behaviors, given by 
Eq. fj4T|) for DRGs, and by an approximation of Eq. (jUJ) for RGs, respectively (see 
the text for details). 
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short time behavior of {0-totif)) in a directed random graph, and in an undirected one 
with the same z = 1.5 average degree. Here, the regime preceding the attainment of 
the asymptotic constant value is relatively long, and one can appreciate the difference 
between the DRG case, described by Eq. (EH), hence (9 to t(t)) ~ z*, and the RG one, 
where at not too short times the data are better in agreement with (6tot{t)) — z t + z t ~ 2 , 
which is an approximation of Eq. (jUJ). Notice moreover that one gets compatible 
asymptotic values in both of the cases. On the right we plot (0 tot (t)), as a function 
of t, in directed random graphs with different z > 1 values, by comparing the curves 
with the expected z l short time behaviors. This figure makes also evident that the 
asymptotic total number of impaired firms is an increasing function of the average 
in/out- degree; moreover, already for z as small as z = 5, it practically coincides with 
the system size N, and the damage spreads in a few steps, t trans ~ log N/ log z, to all 
the firms. 

The results on the average asymptotic density, (p aS ym(z)), as a function of the 
average in/out-degree z, are presented in [Fig. 16], and they turn out to be in 
very good agreement with the expected curve (see Eq. fH3|) ). obtained by evaluating 
numerically the solution v(z) of Eq. (1381) . The data shown in the figure correspond 
to the case of directed random graphs, but we checked that in undirected structures 
with the same average degree one gets definitely compatible results. This agreement 
does also imply that for the considered number of nodes, N = 10000, the corrections 
to the behavior in the limit N — > oo of (p a sym(z)) are practically negligible. 
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Figure 16: Data on free models on directed random graphs, with large network size, 
N = 10000, and deterministically chosen delays all equal to To = 1 day. Here the 
results are averaged over at least J\f s = 200 different configurations of the links. We 
plot the asymptotic average density, (p aS ym(z)) = lim t _ 5 . 00 (p(t)), as a function of the 
mean in/out-degree z, compared with the expected behavior, given by Eq. (1431) . 

To conclude this section, the present results suggest that the study of these simple 
BDE systems on more complex random structures could prove useful also for better 
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understanding the topology of the networks themselves, and in particular for eval- 
uating the size of their connected components, which is difficult to be analytically 
computed in the case of realistic directed networks, where the probability distribu- 
tions of the out and in degree do not factorize [17]. A first step in this direction would 
be to consider (directed) "small- world" networks [39], which in some sense interpo- 
late between the topology of the braid chain and the one of the Erdos-Reny random 
graph. 

It is moreover important to stress that the discussed damage spreading in free BDE 
models is not to be confused with other well known examples of damage spreading, 
such as in particular the spreading of epidemic diseases or of the effects of node 
deletions [40j E]. In fact, on the one hand, here we have assumed that the damage 
propagates to all the nodes in the out-component of the concerned one, though with 
possibly different delays, and, on the other hand, a node attained by the consequences 
of the initial damage is not to be considered removed from the network, since it can 
recover. 

4.3 The free models with random delays 

We now turn onto the study of the free BDE model on a directed random graph, 
when the delays are randomly chosen according to Eq. (Tjfl). The system is described 
by a set of equations analogous to Eq. ( 140~1) : 

Xi{t) = m{t) Y[ Xj{t - nj), i = 1, N. (45) 

where the {t^} variables are independently uniformly distributed in the interval 
[ T min, T max]- Hence the time is in units of r min = 1 day, and we will take as usual 
Tmax = 10 days, limiting moreover the analysis to the case in which the duration of 
the initial damage is r c = r m i n . 

In order to understand the dynamics, we start by noticing that in the first time 
step the signal propagates in average to zjr max other firms. In the next time steps, 
one still expects to observe an exponential damage spreading, for networks above the 
transition point — 1. Moreover, the structure is locally tree-like: this means that, 
at enough short times, the signal propagates along each given branch (and sub-branch 
and so on) of the tree independently; correspondingly, since the propagation paths 
of different time lengths are not concurrent, one expects that the random delays will 
turn out first of all in a global rescaling of the time of a factor l/r av , where r av = (Tjj). 

One obtains: 

(M*)> = —*efY f o r t « logA7logz e// . (46) 

Tmax 

The effective average in/out-degree which appears into this law has to approach z 
when z — > 1 + , hence 

Zeff > Z~ff = z. (47) 
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Figure 17: The total number of impaired firms, 6 to t{t), as a function of time, in the free 
model on a directed random graph with a large network size N = 10000 and different 
z values. The data are averaged over at least J\[ a = 200 different configurations both 
of the links and of the random delays, the last being chosen uniformly distributed 
between T min = 1 day and r max = 10 days. We compare the curves for each z 
value with the expected short time lower limit f~(z,t) = f(z~^,t) and upper limit 
f + (z,t) = f(z^ff,t), where / is given in Eq. ( 1501) . See the text for details. 



Nevertheless, since the signal from a given node propagates to its customers up to 
the time t + r max , one also finds an effective increasing of the average in/out-degree, 
which can be easily overestimated: 

°° / z V z 

**ff < *tff = *Y, [ — ) = 1 - z/T ■ 
\ 'max / "I 'max 

Summarizing, at short times, for z > 1, we get: 

/-(M) = f(z; fp t) < (6 tot (t)) < f(zt fp t) = f + (z,t), (49) 

where 

f(Zeff,t) = -^Z t J f y. (50) 
'max 

We compare in [Fig. 17] the data on (9 to t(t)), for different values of the average 
in/out-degree z > 1, and a large system size N = 10000, with the expected upper 
and lower limits on the short time behaviors, given by Eq. (149)) . The curves do always 
lie between f~(z,t) and f + (z,t) in a large time window; they approach f~(z,t) for 
z —> 1 + , and f + {z) for large z ^> 1; in detail, the data at short times nearly coincide 
with f + (z,t) already for z = 3. 

In the large time limit, since in the giant strongly connected component S gc there 
are loops, one expects that at the end the damage will spread to the whole giant out- 
component O, with the same probability l—v(z) that in the model with deterministic 
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Figure 18: The density of fully active firms p(t) as a function of time in the free 
model on a directed random graph with different average connectivities z > 1. Here 
we consider single typical configurations of the links. We compare the results for 
delays deterministically taken all equal to r = 1 day with the ones for a single typical 
configuration of random delays, chosen uniformly distributed between r m j n = 1 day 
and r max = 10 days. We notice that the densities approach exactly the same values 
in the two cases. 

delays. Correspondingly, we verified numerically that one still observes, for large N 
values, an asymptotic average density of fully active firms (p aS ym) = 2v(z) — v 2 (z), 
in agreement with Eq. ( I43j) . Moreover, for not too small N values, one usually gets 
the same p aS ym for a given network configuration both for deterministic equal delays 
and for randomly chosen nj (see [Fig. 18]). This can be qualitatively understood 
because, if the initially attained firm is in a connected component where there are 
no loops, the economy can recover completely in both of the cases, whereas if it is in 
the giant in-component then the damage spreads to the whole out-component of the 
particular considered network. 

Nevertheless, the asymptotic solutions of the considered free BDE systems, for 
the same directed graph configuration, do not need to be coincident when considering 
equal or random delays. This effect becomes evident for relatively small N, near or 
below the crossover point, hence for z < 1, since it is related to the presence of loops 
in connected components which contain an enough large fraction of the nodes, but 
are not the giant ones. More in detail: 

• when the initially attained firm belongs to a component where there are no 
loops, the asymptotic solutions is given by Xi = 1 Mi in both of the cases; 

• when it belongs to the giant in-component, hence containing an enough large 
number of loops, the asymptotic solutions is given by = Mi G O, and 
Xi = 1 Mi G W \ O, hence the solutions are again the same in both of the cases, 
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Figure 19: The density of fully active firms p(t) as a function of time in a given, 
quite peculiar, random directed network configuration with z = 0.525 and N = 100, 
after a perturbation of node i — 1 of duration r c = 1 day. We compare results for 
deterministically chosen delays, all equal to tq = 1 day (lower curves), to the ones 
for a typical configuration of random delays, uniformly distributed between r m j n = 1 
day and r max = 10 day (upper curves): notice that the asymptotic solutions display 
the same behavior of the density, which oscillates around a value definitely smaller 
than one, though the network is well below the crossover point z ~ 1. See the text 
for details. 
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apart possibly for a small fraction of nodes which approaches zero in the limit 
TV -> oc; 

• in intermediate situations, one usually finds periodic solutions which can be 
largely different. 
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Figure 20: The density of fully active firms p(t) as a function of time, after a pertur- 
bation of node % — 1 of duration r c = 1 day, in some quite peculiar random network 
configurations with z = 0.525, which probability is not negligible already for the rel- 
atively large size N = 100 that we considered. In all of the cases, we compare results 
for deterministically chosen equal delays (lower curve) with the ones for a typical con- 
figuration of random unequal delays (upper curve), as in the previous figure. Notice 
that here the densities approach definitely different asymptotic behaviors in the two 
cases. See the text for details. 



Though a detailed analysis of the finite N behavior of free BDE models on directed 
random graphs is beyond the scope of the present paper, we show in [Fig. 19] and 
[Fig. 20] some examples of peculiar configurations of D(N,p), with N = 100 and 
z = 0.525, well below the crossover point. 

We note that the observed behavior does also depend on the position of the firm 
which is initially destroyed, here taken to be % = 1 in all of the cases. We compare 
the data on the density of fully active firms p(t) for deterministically chosen delays 
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all equal to To = 1 day with the ones for a typical configuration of random unequal 
delays, uniformly distributed between T min = 1 day and r max = 10 days; the duration 
of the initial damage is r c = 1 day in both of the cases. The results confirm that 
one can observe different asymptotic average densities, which implies largely differ- 
ent periodical solutions of the system; moreover, the probability of these disordered 
network configurations is not negligible, since we have observed not compatible p aS ym 
values in ~ 1/10 of the considered samples. 

4.4 The forced models 

We conclude this work by presenting some preliminary results on the forced models, 
as defined by Eq. j2|, on a directed random graphs structure D(N,p). The equation 
for the variable Xi of the system Q can be written as: 

Xi(t) = m(t)xi{t - T ) V Y[ Xj {t - T ) z=l,...,JV. (51) 

when the delays are deterministically taken all equal to To = 1 day, and: 

Xi(t) = m{t) Y[ [xi{t - nj) V Xj (t - T itj )} i = 1, N. (52) 

for random delays. We study, as usual, the damage spreading after the initial de- 
struction of a single firm, that we take for simplicity of duration r c = 1 day. 

As soon as the initial attained firm belongs to a component where there are no 
loops, the economy will recover completely as in the previously considered free models, 
hence the asymptotic solution is Xj = 1 Vi. For large N, this happens in particular 
almost surely below the critical point of the directed random graph zf = 1, and with 
probability v(z) for z > z d , where v(z) is the solution of Eq. ( |38|) . 

Nevertheless, because of the external rescue inputs, one expects that also in the 
region z 3> zf, though the damage spreads almost surely {v (z) <C 1) to the whole 
giant connected out-component, whose size is approximately A" in this limit, the 
average fraction of fully active firms in the asymptotic solution is larger then zero, 
i.e. the economy can partially recover. More in detail, for deterministically taken 
equal delays, since each firm recovers both if all of its suppliers did recover and if 
its activity has been already impaired for a duration t , one expects that, in the 
asymptotic solution, in average the firm is fully active one half of the time. This 
means (p aS ym) = 1 — [1 — v (z)} 2 /2 ~ 1/2, for z ^> 1. In the presence of random 
delays, since there is a larger number of concurrent paths of different time lengths, 
and the activity of a given firm can be impaired for a different duration of time 
depending on the good which is lacking, the asymptotic average density is expected 
to be smaller. 

We present in [Fig. 21] our numerical results on the behavior of p(t), as a function 
of time, for a single, typical, configuration of the directed random graph, of large size 
A^ = 10000, and average in/out-degree z = 7 3> zf. In fact, for deterministically 
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Figure 21: The density of fully active firms p(t) as a function of time, after an initial 
perturbation of node % = 1 of duration r c = 1 day, in the forced model on a directed 
random graph with a large network size N = 10000 and average input/output degree 
z = 7 3> zf . Here we consider a single (typical) configuration of the links, and we 
compare the results for the case of delays taken all equal to To = 1 day with the ones 
for a typical configuration of random unequal delays, chosen uniformly distributed 
between T m i n = 1 day and r max = 10 days. 

chosen equal delays the asymptotic density oscillates around the average value 1/2, 
whereas for random delays it takes the lower average value p aS ym — 0.3. Moreover, in 
the first case the oscillations of the density are quite strong, with lower value ~ 0.4 
and upper value ~ 0.6; hence a finite fraction of the firms, of order 0.2iV ~ 2000, is 
involved into them. 

One can explain qualitatively this result by arguing that the damage starts to 
spread from some firm in the periphery of the connected component of the graph, 
rapidly propagating into the center. Then, because of the external inputs, the most 
of the economy recovers, apart from a few firms which are once again in the periphery, 
and have been reached later by the wave of damage: these will be responsible for the 
forthcoming impulse. Roughly speaking, the nodes with larger in/out-degree are the 
ones occupying the more central positions, whereas a small in-degree means that the 
node is more unlikely reached by the damage propagation and a small out-degree 
means that it is more unlikely to transmit the damage, so that these nodes are to 
be considered in the periphery in both of the cases. This picture is sketched in 
[Fig. 22], and it is quite intriguing, since it implies that there is a large factor of 
unpredictability in the behavior of the solution of the BDE system. A first step 
towards a more quantitative analysis could be to look at the modular structure of the 
directed random graph jH]. Notice however that the random delays have the effect 
of smoothening this oscillating behavior. 
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Figure 22: A qualitative sketch of the expected behavior of the probability distribution 
of impaired firms as a function of the distance of their position from the center of the 
connected component of the graph. See the text for details. 

5 Conclusions 

We have studied the propagation of an initial damage, which destroys a single firm for 
a given time, in production networks. In the considered free models, which represent 
isolated networks, since the local firm dynamics is controlled by a logical AND func- 
tion of the inputs, the damage can invade a finite fraction of the nodes or the whole 
network when these two conditions are fulfilled: i) the nodes mean input connectivity 
is larger than one; ii) the duration of the initial damage is larger than the smallest 
propagation time between two nodes. In particular, for a braid chain (deterministic 
and connected) topology, these conditions mean that in the asymptotic solution the 
production of all the firms is impaired. 

Damage spreading velocity strongly depends upon network topology: we have 
shown that the number of attained firms increases linearly with time for the braid 
chain and exponentially for the directed random graphs. We have also introduced a 
distribution of delays, and we have found that the propagation velocity is dominated 
by the fastest segments (the smallest delay) in the braid chain, while the average 
delay limits speed on the directed random graph. Moreover, in random networks the 
saturation level of the fraction of impaired firms does not depend upon the distribution 
of delays, but only upon the particular network topology through the size of the 
connected component, which is in average not negligible as soon as the mean in/out- 
degree is larger than one. 

External supplies, modeled here in terms of forced networks, limit damage and 
the asymptotic solutions are periodic: waves of damage move across the structure. 
In the considered forced models with distribution of delays, the transient before that 
these solutions are reached diverges exponentially with the network size, though the 
average density of fully active firms approaches a nearly constant value after a defi- 
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nitely shorter effective transient. This effective transient corresponds approximately 
to the duration of the first cycle of propagation of the damage across the connected 
component of the networks, which happens similarly to that in the free models, with 
linear or exponential speed depending on the topology. Finally, periodic dynamics 
are obtained also when the duration of the initial damage is smaller than the smallest 
propagation time between two nearest neighbor nodes. 

The models that we have used are extremely simplified with respect to real pro- 
duction networks, and firms behavior is represented in a very idealized way. Still, 
results from this simple analysis suggest that: 

• Damage spreads in absence of sufficient stocks and flexibility; as a result, pro- 
duction shortages are extended and can survive for periods much longer than 
the duration of the initial local damage. This result suggests that an affected 
economy can suffer from disaster consequences even after all physical damages 
have been repaired. 

• The presence of multiple concurrent production paths with different time lengths 
does not necessarily imply a slowing down of the speed of the signal, which can 
propagate as fast as the shortest path, depending on the topology. This result 
suggests that the most vulnerable supply chains control the macroeconomic 
losses and that vulnerability analysis should focus on the identification of key 
vulnerabilities. 

Appendix A 

The average total number of impaired firms (6tot{t)), in the free model on a braid chain 
network with n — 1, defined accordingly to Eq. flH]) and Eq. fl6]), can be expressed as 
function of the probability V(Tk, k), where Tk is the sum of k random variables r^-: 



TV TV 



<M0> 



i=l i=l J 



k=[t/T m ax]+l i= l k=[t/r max ] + l 

Here [v] means the integer part of v, and we are assuming for simplicity that the 
duration of the initial damage is r c = T min = 1 day, hence we are using the conditions 
on the initial state of the system, {xi(t) = <5i,i(0)} for t G [0, 1]. 
One has, from Eq. (jlj): 



/7~max 



j - 2 



(54) 
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a n, 3 = g 4 , (55) 
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for the mean value and the variance of the distribution of the delays, respectively. For 
k > 2, assuming that the variables are independent (which is clearly an approximation 
for k > N), one can apply the central limit theorem: 



, 1 exp (- SkJ^y ] 

^Tikal 



By using this approach, we find: 



(out)) 



' l t e [0,1) 

V(l, 1) = l/r max te[M) 

i/w + EiS"^ ^g(M, fc) * e [[*], M + 1), 2 <t < r max ^7) 

EE[^:L 1+ i Pg(M, fc) * e [[*], M + i), t> r max 



where the sums are clearly dominated by the terms corresponding to k*(rij) ~ [t], i.e. 
the same that are best approximated by the Gaussian distribution. The corresponding 
average density can be straightaway obtained from Eq. Q). 

In the large t limit, (9 to t(t)) approaches quite rapidly a nearly constant value, that 
in the numerically studied case of r max = 10 is found to be: 

lim (6 tot (t)) ~ lim V V G ([t], k) ~ 0.182, (58) 

k=[[t]/T max ]+l 



which gives, again from Eq. 



lim(p(t)) = l-H^. (59) 



Appendix B 



By using the Boolean rule a A b = a V b, Eq. (II ip . defining the free model on the 
circulant matrix, turns out to be equivalent to: 

n 

%i(t) = y^%i-j(t - T i,i-j)i ( 60 ) 

which, in the case of equal delays {ry = Tq Vi, j}, allows the simplification: 

n n 

^(*) = Yl Yl Xi-h-j 2 (t - 2r ) = 
ii=i j2=i 

ran n 

= ' ' ' Xi -h-32---jk(t ~ ^ r o) — 

ii=ij2=i ife=i 

= Ayr ), (61) 

j=k 
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since the different terms containing the same argument are redundant. Hence, by 
choosing k — [t/r ] = [t], one gets: 

N [t]n 

Otot(t) = 5>*(*) = ~ W) = - 1) + 1 for t G [[t], [t] + 1), (62) 

where we are assuming [t](n — 1) + 1 < N, and we used the conditions on the initial 
state of the system, {xi-j(t) = 5i-j t i for t G [0,r c )}, by considering for simplicity a 
duration of the initial damage r c = r . This analysis confirms that, as soon as n > 1, 
the system is dissipative, and in particular the total number of impaired firms (a 
constant function in each time step) is linearly increasing with time, with slope n—1. 
Correspondingly, for a finite size N, the asymptotically stable solution {xi = OVi}, and 
the zero limit value of the density, are attained after the time tt ra ns — (N — l)/(n — 1), 
in agreement with Eq. ffTT|) . One can moreover analogously show that, for t c ^ tq, the 
density is still linearly decreasing with time, and the asymptotic solution is the same 
for t c > To, whereas it is periodic of period r for r c < r , as described in Section 3.4. 

In order to work out the same approach for randomly distributed delays, we label 
l h = Y^=i3v- One finds: 

n n n / 

ii=ij2=i jfc=i V 

kn / 

= Yl Yl ( f 

j=k {h,l2,-M-h=i V 

Since the {t^} are independent identically distributed random variables which take 
integer values (in unities of r m j„) in the interval [l,r maa; ], it follows that their sums 
Tk = Y%=Q T i-ih,i-ih+i ^ a ^ e integer values in the interval [k, kr max ] with probability 
'PiTk, k); moreover, terms Xi-j with the same argument are redundant. The average 
over the disorder can be therefore computed as: 

(xi(t)) = J HdrV^Xiit) = 

= E E r(rk,k)x M (t-%). (64) 

j=k %=k 

Then, one can choose in particular k values in the different terms such that t — Tk is 
always in the initial interval; in detail, looking for simplicity only at t values which 
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(63) 
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are integer multiple of r max , one has: 



nt 



n(i-l) 

j=t j=t-i 



Tit / Tmax 

... + ^t(t,t/r max ) xi-j(0) 

j—t/Tmax 

t nj 

= E Ew^-^' 

j = t j Tmax — J 

and therefore, up to the time £t rans at which {9 to t(t t rans)) — N, one has: 



(65) 



Hot 



(*)>= E [(n-l)* + W,Ar). 



(66) 



Here, Vt(t,k) is the probability for the signal to have propagated fc steps along the 
chain at time t, which can be approximated by a Gaussian with mean value k(r it j) 
and variance kafj (see Eq. (15B"|) ). but needs to be correctly normalized in order to 
get: 



E v t(t,k) = i 



vt. 



(67) 



In detail, we take Vt{t,k) ~ C(t)VG(t,k), where the normalization constant C(t) 
approaches rapidly its large time limit ~ 1/0.182, in agreement with Eq. ( |58|) . 

In fact, we notice that one would get straightaway a linear behavior of (9tot{t)) 
from Eq. (166]) . by taking Vt(t, k) = 1; nevertheless, the Gaussian approximation and 
the use of the correct normalization constant C(t) are essential in order evaluate 
accurately enough the slope. These results therefore confirm, from a different point 
of view, the analysis of Section 3.5; in particular, for the considered case of n = 20 
and r max = 10, the density (pit)) = 1 — (9 tot (t))/N obtained with this approach is in 
very good agreement with the numerical data, as shown in [Fig. 6]. 



Appendix C 

Here we briefly recall the framework of probability generating functions ^3HH]: this 
allows to simply evaluate the size of connected components in Erdos-Renyi random 
graphs, and it is moreover well suitable to be generalized to directed networks and to 
more complex probability distribution V(k) of the node degrees. 
One defines: 

oo 

£o(y) = E p ^' ( 68 ) 

A;=0 
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which implies: 



(69) 

Analogously, the generating function for the probability distribution Vi(k), to have 
k edges leaving a node apart the one from which the signal arrived, is given by: 



CM VWH„> T.T.oi.k + l)V(k + Dp" ff {y) 

fh(y) - X>(% = — • <*» 

and in the Erdos Renyi random graph one has: 

GM = GM = e zl - y - 1] . (7i) 

One can show in particular that, if Q is the generating function for the probability 
of some property of an object, then the probability of the sum of the same property 
on I independent objects is generated by Q l . Hence, if one defines %o(x) (respectively 
T-Li(x)) as the generating function for the distribution probability of the sizes of the 
connected components (respectively of the clusters), i.e. of the nodes which can be 
reached from a randomly chosen one (respectively from the end of one of the edges 
of a randomly chosen one), it is: 



H^y) = yX)7>i(*)[Hi(v)]* = vGi[H 1 { 

k=0 

oo 

-H (y) = ^2V(k)[HM k =yg [HM. (72) 

k=0 

Therefore, one obtains a set of equations which can be in principle solved self- 
consistently. Most importantly from our present point of view, one gets immediately 
the mean connected component size: 

• below the transition one can show that 

(s) = H'o(l) = 1 + G'o(m' 1 (l) = 1 + - G '°}P„„ (73) 



G'i(l) 



which means, in the Erdos-Renyi random graph, 



as expected. 



above the transition, T-Lq{x) and 1-Li{x) can be defined as the generating function 
for the distribution probability of the finite size connected components, which 
have still a tree- like structure; correspondingly H (l) = v < 1, and one has 

v = i-v = i-g (v) 

v = g x {y) 

which for Erdos-Renyi random graphs gives v = e z( - v ~ l \ i.e. Eq. 
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